clear
* Change directory to where the folder with the data is
cd "Data Archive TWT Predictions\KTHDP"

*Load the dataset
use DatasetKTHDP.dta,clear

*Summarize the dataset
sum prob1 prob2 out1 out2 choiceleft
tab subj
tab choiceleft
*transform choice dummy from 0,1 to -1,1
replace choiceleft=-1 if choiceleft==0
tab choiceleft
save BaseEst, replace



*Start the inizialization of the ML program
clear
program drop _all
program define ML_eut0

	* specify the arguments of this program
	args lnf r noise
		
	* declare the temporary variables to be used
	tempvar euL euR uL1 uL2 uR1 uR2 choicea probl prob2 out1 out2 maxOut euDiff
	
	* please do not display all these steps on the screen for every iteration
	quietly {

	
	*choiceA probAl probBl outA1 outB1 maxOut probA2 probB2 outA2 outB2
		generate int    `choicea' = $ML_y1
		generate double `probl' = $ML_y2 
		generate double `prob2' = $ML_y3 
		generate double `out1' = $ML_y4
		generate double `out2' = $ML_y5
		generate double `maxOut' = $ML_y6

		*calculate utility each option usinf CARA, plug in here different utility functions (e.g. CRRA) for alternative estimations. 
		*for RPM one just needs to ocmpare the risk attitude which would make a DM indifferent for each lottery pair with the current best estimate of the risk attitude. Can skip the next two passages.
		generate double `uL1' = ((1-exp(-`r'*`out1'))/(1-exp(-`r'*(`maxOut'))))
		generate double `uR1' = ((1-exp(-`r'*`out2'))/(1-exp(-`r'*(`maxOut'))))
		
		*calculate EU of each lottery		
		generate double `euL' = (`probl'*`uL1')
		generate double `euR' = (`prob2'*`uR1')

		*difference EU divided by noise
		generate double `euDiff' = (`euL' - `euR')/`noise'

		*evaluate the likelihood choice based on EU diff
		replace `lnf' = ln(normal( `euDiff')) if `choicea'==1
		replace `lnf' = ln(normal(-`euDiff')) if `choicea'==-1
	}

	
end



*Loop estimation over participants, skip those subjects for whom the estimation do not converge.
local i = 1
while `i' < 31 {
local j = 1
while `j' < 6 {
display `i'

use baseest, clear
keep if subj==`i'
* This is an out-of-sample estimation, risk parameters are estimated from a subset of choices to predict a different set. 
keep if idleft==`j' | idright==`j'

egen m1 = max(out1)
egen m2 = max(out2)
gen maxOut=max(m1,m2)
drop m1 m2
 
gen choicea=choiceleft
sum choicea prob1 prob2 out1 out2 maxOut
gen est_r_u`j'=.
gen noise_u`j'=.

display `i'
capture{
ml model lf ML_eut0 (r: choicea prob1 prob2 out1 out2 maxOut  = ) (noise: ), display difficult  maximize
matrix b = e(b)  
matrix list b  
replace est_r= b[1,1]  
replace noise = b[1,2]
matrix drop b
  }
*keep only subjects' identifier, risk attitudes, and noise parameter.
keep  subj est_r_u*  noise_u*
collapse subj est_r_u`j'  noise_u`j'
save individualiskCARAsubj`i'u`j', replace
local j = `j' + 1
}
local i = `i' + 1
}

*Merge individual datasets into one.
local i = 1
while `i' < 31 {
use  individualiskCARAsubj`i'u1, clear
local j = 2
while `j' < 6 {
append using  individualiskCARAsubj`i'u`j'
local j = `j' + 1
}
save individualiskCARAsubj`i', replace
local i = `i' + 1
}
clear 
local i = 1
while `i' < 31 {
append using individualiskCARAsubj`i'
local i = `i' + 1
}
collapse est_r_u*  noise_u*, by(subj)
save UtilityEstimationCARAKTHDP, replace
*Repeat the same procedure for RUM-CRRA, RPM-CARA, and RPM-CRRA.
*saved files should be:
*	UtilityEstimationCRRAKTHDP
*	UtilityEstimationCARARPMKTHDP
*	UtilityEstimationCRRARPMKTHDP
